#-----------------------------------#
# Replication of:
# Seeking Opportunity in the Knowledge Economy: Moving Places, Moving Politics?
# Valentina Consiglio, Thomas Kurer

# File: 7a_models_mechanisms

# Note: This file codes the models included in the mechanisms section of the paper. It takes the following steps:
# (1) Analysis of symmetry of effects (up- vs. downward moves)
# (2) Analysis of economic and cultural mechanisms
# (3) Analysis of symmetry of economic and cultural mechanisms
# The data is loaded and prepared in 5_models_main_NUTS3; ensure the script was run before this file. 
#-----------------------------------#


#-------------------------------------------------------------------------------# 
#### Preparation ####
#-------------------------------------------------------------------------------# 

## Texreg input 
vec_gof_names <- c("Num. obs." = "N",
                   "Num. groups: pid" = "N individuals",
                   "Num. groups: syear" = "N years")


#-----------------------------------------------------------------------------------------------------#
#### 1 | Symmetry of Effects ####
#-----------------------------------------------------------------------------------------------------#

#-----------------------------------------------------------------------------------------------------#
##### 1.1 | Upward move  #####
#-----------------------------------------------------------------------------------------------------#

#-----------------------------------------------------------------------------------------------------#

## Opportunity moves and political integration and orientation ##

# Without multiple movers


# Unbalanced
fe_main_polint_polor2_upmove_ub2_only <-  list(
  "Volunteering (Yes/No)" = feols(ft_volunt_d ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  
  "Local political activism (Yes/No)" = feols(ft_localpol_d ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  
  "Voting in federal election (Yes/No)" = feols(votefedel_lead ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009_upmove_ub2),
  
  "Party leaning (Yes/No)" = feols(partyleand2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| 
                                     pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  
  "Party leaning intensity (0-5)" = feols(partyleanint ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| 
                                            pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  
  "Political orientation (0 L – 10 R)" = feols(polor ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                 pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009_upmove_ub2),
  
  "Centre left party ID" = feols(centre_left2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                   pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  
  "Centre right party ID" = feols(centre_right2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2)
)


texreg(fe_main_polint_polor2_upmove_ub2_only)

texreg(fe_main_polint_polor2_upmove_ub2_only, digits = 3,
       label = "tab:fe_main_polint_polor2_upmove_ub2_only",
       caption = "Upward opportunity moves and political integration and orientation",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. All models are restricted to one-time upward movers and include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E30.tex"))



fe_main_partid2_upmove_ub2_only <-  list(
  "SPD" = feols(spd_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  "CDU/CSU" = feols(cducsu_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                      pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  "Gruene" = feols(greens_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                     pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  "FDP" = feols(fdp_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  "Linke" = feols(linke_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_upmove_ub2),
  "AfD" = feols(afd_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2014_upmove_ub2)
)

texreg(fe_main_partid2_upmove_ub2_only)


texreg(fe_main_partid2_upmove_ub2_only, digits = 3,
       label = "tab:fe_main_partid2_upmove_ub2_only",
       caption = "Upward opportunity moves and party identification",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                            
       custom.note = paste("%stars. All models are restricted to one-time upward movers and include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E31.tex"))



#-----------------------------------------------------------------------------------------------------#
##### 1.2 | Downward move  #####
#-----------------------------------------------------------------------------------------------------#

# Unbalanced
fe_main_polint_polor2_downmove_ub2_only <-  list(
  "Volunteering (Yes/No)" = feols(ft_volunt_d ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  
  "Local political activism (Yes/No)" = feols(ft_localpol_d ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  
  "Voting in federal election (Yes/No)" = feols(votefedel_lead ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009_downmove_ub2),
  
  "Party leaning (Yes/No)" = feols(partyleand2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| 
                                     pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  
  "Party leaning intensity (0-5)" = feols(partyleanint ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| 
                                            pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  
  "Political orientation (0 L – 10 R)" = feols(polor ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                                 pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009_downmove_ub2),
  
  "Centre left party ID" = feols(centre_left2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                   pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  
  "Centre right party ID" = feols(centre_right2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2)
)


texreg(fe_main_polint_polor2_downmove_ub2_only)

texreg(fe_main_polint_polor2_downmove_ub2_only, digits = 3,
       label = "tab:fe_main_polint_polor2_downmove_ub2_only",
       caption = "Downward opportunity moves and political integration and orientation",
       custom.coef.map = list('opportunity_pc01_wn_rc' = 'Opportunity Index (rev)'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. All models are restricted to one-time downward movers and include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E33.tex")
)



fe_main_partid2_downmove_ub2_only <-  list(
  "SPD" = feols(spd_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  "CDU/CSU" = feols(cducsu_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                      pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  "Gruene" = feols(greens_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                     pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  "FDP" = feols(fdp_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  "Linke" = feols(linke_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_downmove_ub2),
  "AfD" = feols(afd_co2 ~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                  pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2014_downmove_ub2)
)

texreg(fe_main_partid2_downmove_ub2_only)


texreg(fe_main_partid2_downmove_ub2_only, digits = 3,
       label = "tab:fe_main_partid2_downmove_ub2_only",
       caption = "Downward opportunity moves and party identification",
       custom.coef.map = list('opportunity_pc01_wn_rc' = 'Opportunity Index (rev)'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                            
       custom.note = paste("%stars. All models are restricted to one-time downward movers and include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E34.tex"))



#-------------------------------------------------------------------------------# 
#### 2 | Economic and Cultural mechanisms ####

#-------------------------------------------------------------------------------# 

##### 2.1 | Economic #####


# Main #
econ_DVs <- c("ln_pglabgro2", "pgisei08", "worr_wpsec")

names_econ_DVs <- c(ln_pglabgro2 = "Log gross monthly earnings", 
                    pgisei08 = "ISEI-score",
                    worr_wpsec = "Worries about workplace security"#, 
                    #income = "Equiv. HH income"
)


fe_mech_econ <- list()

# Loop over econ CVs
for (econ_DV in econ_DVs) {
  
  # Define model
  formula <- as.formula(paste(econ_DV, "~ opportunity_pc01_wn + pgexpft + pgexpft^2 + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_econ[[names_econ_DVs[[econ_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "krs", data = df_2009)
}

texreg(fe_mech_econ)
texreg(fe_mech_econ, digits = 3,
       label = "tab:fe_mech_econ",
       caption = "Opportunity moves and socio-economic outcomes",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_econ)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_econ))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. Models with socio-economic outcomes include work experience, work experience squared, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Table_5a.tex"))

# Moved only

fe_mech_econ_moved <- list()

# Loop over econ CVs
for (econ_DV in econ_DVs) {
  
  # Define model
  formula <- as.formula(paste(econ_DV, "~ count_moved_krs2 + pgexpft + pgexpft^2  + as.factor(educ_h) + as.factor(hhtype) | pid + syear"))
  # + as.factor(educ_h) + as.factor(hhtype)
  # Run regression 
  fe_mech_econ_moved[[names_econ_DVs[[econ_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "krs", data = df_2009_womm)
}

texreg(fe_mech_econ_moved)

texreg(fe_mech_econ_moved, digits = 3,
       label = "tab:fe_mech_econ_moved_womm",
       caption = "Kreis-level moves: Socio-economic Outcomes",
       custom.coef.map = list('count_moved_krs2' = 'Kreis-level move'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_econ_moved)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_econ_moved))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. Kreis-level move is defined as a cumulative count-variable. All models include relevant control variables: Models with socio-economic outcomes include work experience, work experience squared, education group, and household type as control variables. Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E36a.tex"))



#------------------------------------------------------------------------------#
##### 2.2 | Cultural #####
#------------------------------------------------------------------------------#

# DVs on cultural variables

cult_DVs <- c("ft_opera_d", "ft_cinema_d","ft_artmus_d") # 

names_cult_DVs <- c(ft_opera_d = "Cultural events (dummy)", 
                    ft_cinema_d = "Cinema, concerts, clubs (dummy)", 
                    ft_artmus_d = "Artistic/musical activities (dummy)"
)

fe_mech_culture <- list()

# Loop over cult CVs
for (cult_DV in cult_DVs) {
  
  # Define model
  formula <- as.formula(paste(cult_DV, "~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_culture[[names_cult_DVs[[cult_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "krs", data = df_2009)
}


texreg(fe_mech_culture)
texreg(fe_mech_culture, digits = 3,
       label = "tab:fe_mech_culture",
       caption = "Kreis-level moves: Cultural Outcomes",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_culture)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_culture))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. Models with cultural activities outcomes include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Table_5b.tex"))


# moved (stay treated)

fe_mech_culture_moved <- list()

# Loop over cult CVs
for (cult_DV in cult_DVs) {
  
  # Define model
  formula <- as.formula(paste(cult_DV, "~ count_moved_krs2 + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_culture_moved[[names_cult_DVs[[cult_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "krs", data = df_2009_womm)
}


texreg(fe_mech_culture_moved)
texreg(fe_mech_culture_moved, digits = 3,
       label = "tab:fe_mech_culture_moved",
       caption = "Kreis-level moves and cultural activities",
       custom.coef.map = list('count_moved_krs2' = 'Kreis-level move'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_culture_moved)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_culture_moved))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. Kreis-level move is defined as a cumulative count-variable. All models include relevant control variables: Models with cultural activity outcomes include age group, education group, and household type as control variables. Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E36b.tex"))

#-----------------------------------------------------------------------------#
#### 3 | Symmetry mechanisms ####
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
###### 3.1 | Upward Econ ######
#------------------------------------------------------------------------------#


# Main #
econ_DVs <- c("ln_pglabgro2", "pgisei08", "worr_wpsec")

names_econ_DVs <- c(ln_pglabgro2 = "Log gross monthly earnings", 
                    pgisei08 = "ISEI-score",
                    worr_wpsec = "Worries about workplace security"
                    )

fe_mech_econ_upmove_ub2 <- list()

# Loop over econ CVs
for (econ_DV in econ_DVs) {
  
  # Define model
  formula <- as.formula(paste(econ_DV, "~ opportunity_pc01_wn + pgexpft + pgexpft^2 + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_econ_upmove_ub2[[names_econ_DVs[[econ_DV]]]] <- feols(formula, panel.id = ~pid + syear, 
                                                                cluster = "krs", data = df_2009_upmove_ub2)
}


texreg(fe_mech_econ_upmove_ub2)


texreg(fe_mech_econ_upmove_ub2, digits = 3,
       label = "tab:fe_mech_econ_upmove_ub2",
       caption = "Upward opportunity moves: Socio-economic Outcomes",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_econ_upmove_ub2)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_econ_upmove_ub2))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. The sample is restricted to upward moves only. All models include work experience, work experience squared, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E32a.tex"))


#------------------------------------------------------------------------------#
###### 3.2 | Downward Econ ######
#------------------------------------------------------------------------------#

# Main #
econ_DVs <- c("ln_pglabgro2", "pgisei08", "worr_wpsec")

names_econ_DVs <- c(ln_pglabgro2 = "Log gross monthly earnings", 
                    pgisei08 = "ISEI-score",
                    worr_wpsec = "Worries about workplace security"#, 
                    #income = "Equiv. HH income"
)


fe_mech_econ_downmove_ub2 <- list()

# Loop over econ CVs
for (econ_DV in econ_DVs) {
  
  # Define model
  formula <- as.formula(paste(econ_DV, "~ opportunity_pc01_wn_rc + pgexpft + pgexpft^2 + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_econ_downmove_ub2[[names_econ_DVs[[econ_DV]]]] <- feols(formula, panel.id = ~pid + syear, 
                                                                  cluster = "krs", data = df_2009_downmove_ub2)
}


texreg(fe_mech_econ_downmove_ub2, digits = 3,
       label = "tab:fe_mech_econ_downmove_ub2",
       caption = "Downward opportunity moves: Socio-economic Outcomes",
       custom.coef.map = list('opportunity_pc01_wn_rc' = 'Opportunity Index (rev)'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_econ_downmove_ub2)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_econ_downmove_ub2))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. The sample is restricted to downward moves only. All models include work experience, work experience squared, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2010-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E35a.tex"))

#------------------------------------------------------------------------------#
###### 3.3 | Upward Cult ######
#------------------------------------------------------------------------------#


cult_DVs <- c("ft_opera_d", "ft_cinema_d","ft_artmus_d")

names_cult_DVs <- c(ft_opera_d = "Cultural events (dummy)", 
                    ft_cinema_d = "Cinema, concerts, clubs (dummy)", 
                    ft_artmus_d = "Artistic/musical activities (dummy)"
)

fe_mech_culture_upmove_ub2 <- list()

# Loop over cult CVs
for (cult_DV in cult_DVs) {
  
  # Define model
  formula <- as.formula(paste(cult_DV, "~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_culture_upmove_ub2[[names_cult_DVs[[cult_DV]]]] <- feols(formula, panel.id = ~pid + syear, 
                                                                   cluster = "krs", data = df_2009_upmove_ub2)
}


texreg(fe_mech_culture_upmove_ub2, digits = 3,
       label = "tab:fe_mech_culture_upmove_ub2",
       caption = "Upward opportunity moves: Cultural Outcomes",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_culture_upmove_ub2)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_culture_upmove_ub2))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. The sample is restricted to upward moves only. All include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E32b.tex"))



#------------------------------------------------------------------------------#
###### 3.4 | Downward Cult ######
#------------------------------------------------------------------------------#

# Dummies
cult_DVs <- c("ft_opera_d", "ft_cinema_d","ft_artmus_d")

names_cult_DVs <- c(ft_opera_d = "Cultural events (dummy)", 
                    ft_cinema_d = "Cinema, concerts, clubs (dummy)", 
                    ft_artmus_d = "Artistic/musical activities (dummy)"
)


fe_mech_culture_downmove_ub2 <- list()

# Loop over cult CVs
for (cult_DV in cult_DVs) {
  
  # Define model
  formula <- as.formula(paste(cult_DV, "~ opportunity_pc01_wn_rc + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_culture_downmove_ub2[[names_cult_DVs[[cult_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "krs", data = df_2009_downmove_ub2)
}




texreg(fe_mech_culture_downmove_ub2, digits = 3,
       label = "tab:fe_mech_culture_downmove_ub2",
       caption = "Downward opportunity moves: Cultural Outcomes",
       custom.coef.map = list('opportunity_pc01_wn_rc' = 'Opportunity Index (rev)'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_culture_downmove_ub2)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_culture_downmove_ub2))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. The sample is restricted to downward moves only. All models include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E35b.tex"))



#-------------------------------------------------------------------------------# 
#### 4 | Occupational stayers ####

#-------------------------------------------------------------------------------# 

# Check variables

table(df$occ_chg, useNA =  "always")
table(df$occ_chg_move, useNA =  "always")


pid_with_na <- df %>%
  filter(is.na(occ_chg_move)) %>%
  pull(pid)


# Those NA who have NA on isco1d --> code them as 0. 
df_2009$occ_chg_move[is.na(df_2009$occ_chg_move)] <- 0
df$occ_chg_move[is.na(df$occ_chg_move)] <- 0
df_2014$occ_chg_move[is.na(df_2014$occ_chg_move)] <- 0

table(df$isco_h_1d, useNA =  "always")

# Check number/share of moves that are also occupational moves
table(df$occ_chg_move, df$moved_krs2, useNA =  "always")

df_sample <- df %>%
  filter(!is.na(spd_co2),
         !is.na(opportunity_pc01_wn),
         !is.na(nuts1),
         !is.na(age_grp3),
         !is.na(educ_h),
         !is.na(hhtype)
  )

# Check occupational movers / stayers
table(df_sample$occ_chg_move, df_sample$moved_krs2, useNA =  "always")

## Subset to occupational changers ##
 # i.e. drop if move but doesn't change

df <- df %>%
  group_by(pid) %>%
  mutate(flag_occstay = ifelse(any(moved_krs2 == 1 & occ_chg_move == 0), 1, 0),
         flag_occmove = ifelse(any(moved_krs2 == 1 & occ_chg_move == 1), 1, 0)
  ) %>%
  ungroup()

table(df$flag_occstay, df$flag_occmove, useNA = "always")
table(df$flag_occstay, useNA = "always")
table(df$flag_occmove, useNA = "always")

df_2009 <- df_2009 %>%
  group_by(pid) %>%
  mutate(flag_occstay = ifelse(any(moved_krs2 == 1 & occ_chg_move == 0), 1, 0),
         flag_occmove = ifelse(any(moved_krs2 == 1 & occ_chg_move == 1), 1, 0)
  ) %>%
  ungroup()

df_2014 <- df_2014 %>%
  group_by(pid) %>%
  mutate(flag_occstay = ifelse(any(moved_krs2 == 1 & occ_chg_move == 0), 1, 0),
         flag_occmove = ifelse(any(moved_krs2 == 1 & occ_chg_move == 1), 1, 0)
  ) %>%
  ungroup()


#-------------------------------------------------------------------------------# 
## Political integration and orientation ##
#-------------------------------------------------------------------------------# 

fe_main_polint_polor2_subset_occstay <-  list(
  "Volunteering (Yes/No)" = feols(ft_volunt_d ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                  | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  
  "Local political activism (Yes/No)" = feols(ft_localpol_d ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                              | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  
  "Voting in federal election (Yes/No)" = feols(votefedel_lead ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                                | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009 %>% filter(flag_occmove != 1)),
  
  "Party leaning (Yes/No)" = feols(partyleand2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                   | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  
  "Party leaning intensity (0-5)" = feols(partyleanint ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                          | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  
  "Political orientation (0 L – 10 R)" = feols(polor ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                               | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2009 %>% filter(flag_occmove != 1)),
  
  "Centre left party ID" = feols(centre_left2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                 | pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  
  "Centre right party ID" = feols(centre_right2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) |
                                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1))
)


texreg(fe_main_polint_polor2_subset_occstay)


texreg(fe_main_polint_polor2_subset_occstay, digits = 3,
       label = "tab:fe_main_polint_polor2_occstay",
       caption = "Opportunity moves and political integration and orientation (occupational stayers)",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. The treatment is restricted to those who relocate and do not change their 1-digit ISCO occupation. All models include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E37.tex"))


#--------------------------------------------------------------------------------------------#

#-------------------------------------------------------------------------------# 
## Party identification ##
#-------------------------------------------------------------------------------# 

# - alternative coding where all those with 0 party identification are in the reference category of the dummies

fe_main_partid2_subset_occstay <-  list(
  "SPD " = feols(spd_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  "CDU/CSU " = feols(cducsu_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                       pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  "Gruene " = feols(greens_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                      pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  "FDP " = feols(fdp_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  "Linke" = feols(linke_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                    pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df %>% filter(flag_occmove != 1)),
  "AfD " = feols(afd_co2 ~ opportunity_pc01_wn + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "krs", data = df_2014 %>% filter(flag_occmove != 1))
)

texreg(fe_main_partid2_subset_occstay)


texreg(fe_main_partid2_subset_occstay, digits = 3,
       label = "tab:fe_main_partid2_subset_occstay",
       caption = "Opportunity moves and party identification (occupational stayers)",
       custom.coef.map = list('opportunity_pc01_wn' = 'Opportunity Index'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                            
       custom.note = paste("%stars. The treatment is restricted to those who relocate and do not change their 1-digit ISCO occupation. Party leanings include coalitions in the political direction of the respective party (i.e., including left coalitions for SPD, Gruene, and Linke; right coalitions for CDU/CSU, FDP, and AfD). All models include age group, education group, and household type as control variables; Standard errors are clustered at the Kreis-level. Source: SOEP v.37, 2010-2020."),
       file = paste0(tables_path, "mechanisms/Appendix_Table_E38.tex"))

#--------------------------------------------------------------------------------------------#


#-------------------------------------------------------------------------------# 

# Remove all model lists from storage

rm(list = ls(pattern = "^fe_"))
#-----------------------------------------------------------------------------------------------------#


